Structural relaxation in supercooled orthoterphenyl 
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We report molecular-dynamics simulation results performed for a model of molecular liquid or- 
thoterphenyl in supercooled states, which we then compare with both experimental data and mode- 
coupling-theory (MCT) predictions, aiming at a better understanding of structural relaxation in 
orthoterphenyl. We pay special attention to the wavenumber dependence of the collective dynam- 
ics. It is shown that the simulation results for the model share many features with experimental 
data for real system, and that MCT captures the simulation results at the semiquantitative level 
except for intermediate wavenumbers connected to the overall size of the molecule. Theoretical 
results at the intermediate wavenumber region are found to be improved by taking into account 
the spatial correlation of the molecule's geometrical center. This supports the idea that unusual 
dynamical properties at the intermediate wavenumbers, reported previously in simulation studies for 
the model and discernible in coherent neutron-scattering experimental data, are basically due to the 
coupling of the rotational motion to the geometrical-center dynamics. However, there still remain 
qualitative as well as quantitative discrepancies between theoretical prediction and corresponding 
simulation results at the intermediate wavenumbers, which call for further theoretical investigation. 

PACS numbers: 61.20. Ja,61.20.Lc,61.25.Em,64.70.Pf 



o 



> 

O 

o 



I 

O 

o 



X 
J3 



I. INTRODUCTION 

Describing the microscopic origin of structural slow- 
ing down on cooling or compressing glass-forming liquids 
is one of the most challenging problems in condensed 
matter physics. During the past decade the research in 
this field was strongly influenced by the idealized mode- 
coupling theory (MCT) for the evolution of structural 
relaxation [1-3]. The theory predicts the structural ar- 
rest - also referred to as the idealized liquid-glass tran- 
sition - driven by the mutual blocking of a particle and 
its neighbors ( "cage effect" ) at a critical temperature Tc 
which is located above the glass-transition temperature 
Tg. For temperatures close to but above Tc, MCT pre- 
dicts universal scaling laws and power laws for describ- 
ing the glassy slow structural relaxation. Although such 
complete structural arrest at Tc is not observed in ex- 
periments and therefore the idealized theory cannot lit- 
erally be applied for describing dynamics below Tc, ex- 
tensive tests of the theoretical predictions carried out so 
far against experimental and computer-simulation results 
suggest that MCT deals properly with some essential fea- 
tures of supercooled liquids [2, 3]. 

The molecular van der Waals liquid orthoterphenyl 
(OTP) has long been used as a model system in the study 
of the glass transition. (See, e.g., Refs. [4-6] and papers 
quoted therein.) Extensive experiments on OTP have 
been performed to monitor the onset of glassy structural 
relaxation on microscopic time and length scales. Using 
in particular quasielastic neutron scattering, the decay of 
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collective and self density fluctuations has been measured 
as a function of temperature, pressure, and wavenum- 
ber [4, 7-12]. Based on these studies, the validity of the 
universal predictions of MCT, such as the factorization 
theorem and the time-temperature superposition princi- 
ple, has been established. However, there are nonuniver- 
sal aspects in the glassy structural relaxation which can- 
not be elucidated solely by those universal predictions. 
For example, parameters describing the decay of the col- 
lective density fluctuations in the a-relaxation regime ex- 
hibit a characteristic wavenumber dependence, which is 
an important nonuniversal aspect to be accounted for by 
the theory [4, 12]. As will be discussed below, we found 
from the analysis of molecular-dynamics (MD) simula- 
tion results interesting dynamical properties of OTP at 
intermediate wavenumbers, and their theoretical investi- 
gation is one of the main points of the present paper. 

One of the distinctive features of MCT is that it can 
also make predictions concerning nonuniversal aspects, 
such as the value of Tc and the details of the time and 
wavenumber dependence of various dynamical quantities, 
provided the system's static structure factor is known 
with sufficient accuracy. Utilizing this feature, there 
have been quantitative tests of the theory concerning 
nonuniversal as well as universal aspects using as input 
only the static structure factor determined from integral- 
equation theories or from computer simulations [13-19]. 
In the present paper, we apply MCT to discuss in de- 
tail properties of the slow structural relaxation in OTP, 
paying special attention to the wavenumber dependence 
of the structural a relaxation of the collective dynam- 
ics. This will done by regarding MD simulation results 
as a bridge connecting experimental data and theoretical 
predictions. 

Although OTP is one of the simplest molecular sys- 
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terns from an experimental side, it is rather a compli- 
cated system for a theoretical treatment. Therefore, it is 
unavoidable to deal with a simple model for OTP which 
is still efRcient in mimicking the complexity of the dy- 
namical behavior of real system. In this respect, Lewis 
and Wahnstrom (LW) [20] introduced a particularly use- 
ful three-site model, each site playing the role of a whole 
phenyl ring, and this model will be considered in this 
paper. 

For the LW OTP model, static and dynamic properties 
in supercooled states have been extensively studied based 
on MD simulations [21-23]. In the present work, the sim- 
ulation results for the static structure factors serve as in- 
put to the theory, and those for dynamics can be used 
in testing the so-obtained theoretical results. Thereby, 
a stringent test of MCT predictions against the simula- 
tion results can be performed. Furthermore, by making 
connections between the simulation results for LW OTP 
and experimental ones for real system, the relevance of 
our theoretical results and their interpretation in under- 
standing experimental data can be established. 

A theoretical analysis for LW OTP has already been 
presented in Ref. [21] based on a simplified theory which 
is essentially the same as the one for spherical particles. 
But, since OTP is a molecular system, a full molecu- 
lar approach is desirable. Recently, MCT for spherical 
particles has been extended to a theory for molecules. 
This has been done based on the tensor- and site-density 
formulations. In the tensor-density formulation, the 
density-fluctuation correlator is generalized to the one 
of infinite matrices of correlation functions formed with 
tensor-density fluctuations [17-19, 24-29]. However, the 
tensor-density formulation has the difliculty that the re- 
sulting equations are so involved, and it is not obvious 
whether those equations can be numerically solved within 
the regime of glassy dynamics. To overcome this diffi- 
culty, it has been suggested to base MCT for molecular 
systems on the site representation [30-34]. The fluctua- 
tions of the interaction-site densities have been used as 
the basic variables to describe the structure of the sys- 
tem. As a result, the known scalar MCT equations for 
the density fluctuations in simple systems have been gen- 
eralized to n-hy-n matrix equations for the interaction- 
site-density fluctuations, where n denotes the number 
of atoms (or sites) forming the molecule. Thus, rela- 
tively simple equations of motions can be obtained within 
the site-representation, and these MCT equations will be 
solved in the present work to discuss the structural re- 
laxation in OTP. 

The paper is organized as follows. In Sec. II, the LW 
OTP model shall be introduced, and static as well as 
dynamic quantities to be used in discussions of simu- 
lation and theoretical results are defined. In Sec. Ill, 
MD-simulation results for LW OTP arc summarized, and 
possible connections to the experimental results for real 
system are established. In Sec. IV, after reviewing rel- 
evant MCT equations based on the site representation, 
theoretical results are presented and compared with the 
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FIG. 1: Schematic representation of the three-site LW OTP 
molecule in tlie body-fixed molecular XY plane where the 
origin is taken to be the geometrical center (GC), f'p'^ = 
(1/3) ^a,=i ^^'^ ^ along the symmetry axis of 

the molecule. The diameter of each site is taken from the LJ 
parameter a = 0.483 nm. The dotted circle is drawn with the 
van der Waals radius rw = 0.37 nm [37] and taking GC as 
the origin. 

simulation results. The paper is summarized in Sec. V 
with some concluding remarks. 

II. MODEL 

The LW OTP molecule designed by Lewis and Wahn- 
strom [20] is a rigid isosceles triangle; the length of the 
two short sides of the triangle is £ = 0.483 nm and the 
angle between them is ^ = 75°. Each of the three sites 
represents an entire phenyl ring of mass m » 78 amu, 
and is described by a Lennard-Jones (LJ) sphere whose 
interaction potential is given by 

Vir) ^ 4e [(o-/r)i2 _ (^ajrfX + A + Br., (1) 

with e = 5.276 kJ/mol, a = 0.483 nm, A = 0.461 kJ/mol, 
and B = —0.313 kJ/(mol nm). Schematic representation 
of the LW OTP molecule is presented in Fig. 1. The 
shape of the molecule and the LJ parameters have been 
chosen to reproduce some bulk properties of the OTP 
molecule. The values of A and B are selected in such 
a way that the potential and its first derivative are zero 
at the cutoff Tc = 1.2616 nm adopted in MD simula- 
tions. MD simulation results on the thermodynamic and 
dynamic properties of LW OTP have been reported in 
Rcfs. [20-23]. 

Discussions on the dynamics shall be done based on 
the site-site density correlators 

Ffit) = iPWpUO))/N (a, 6 = 1,2, 3), (2) 
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and their derivatives to be defined below. Here p^{t) = 

SjLi6xp[ig • rj'{t)], with r?(t) being the position vec- 
tor of site a in jth molecule at time t, denotes the site- 
density fluctuations, N the number of molecules in the 
system, and (• • • ) canonical averaging for temperature T. 
(As depicted in Fig. 1, we shall use the convention that 
a = 1 refers to the central site and a = 2,3 to the ad- 
jacent sites.) Because of isotropy, F^^{t) depends only 



on the wavenumber q 



The initial values consti- 



tute the site-site static structure factors, S^' 



which provide the simplest information on the equilib- 
rium structure of the system. 

We also introduce tagged-molecule corre- 
lators, F-'M = (pf,,(t)>^-,,(0)), in which 



exp[iq ■ fg{t)] with f^{t) denoting the 
position vector of site a in the tagged molecule (labeled 
s) at time t. The initial value shall be denoted as 
w^^ = F^^(O)- For a rigid molecule, it is given by 
wf = (5"'' + (1 - (5"'')jo(gZ°^), where jo(x) denotes the 
Oth-order spherical Bessel function and 1°''' the distance 
between sites a and b. 

Let us consider a more convenient representation of 
the site-site density correlators which exploits the C2v 
symmetry of the LW OTP molecule. For this purpose, 
we introduce the following density correlators 

F^it) = {P^ityp^m/N for X = N, Z, and Q, (3) 

defined in terms of the linear combinations 

= (Pf + pI + P| = - P} - 4)/^^' 

= (4) 

of the site-density fiuctuations p| (a = 1, 2, or 3). One 
can easily show that F?'{t) can be expressed in terms of 



F;^\t): one finds, for example, F^{t) = Ea,fe=i 
Their normalized correlators shall be denoted as (f)^{t) = 
F^{t)/S^ with the corresponding static structure fac- 
tor = F^{0). We also introduce self-parts of these 
correlators, to be denoted as F^g{t), which are defined 
similarly in terms of the site-density fluctuations pi ^ of 
the tagged molecule. Normalized tagged-molecule's cor- 
relators shall be written as (j^fsit) = F^g{t)/w^ with 
w^ = FlM- 

Due to the C^v symmetry of the LW OTP 
molecule, one finds that density correlators involving pS, 

{p^it)* p^) /N , become nonzero only for X = Q, and the 
nonzero correlator for X = Q is identical to its self-part, 
i.e., F^(t) = F^g{t). Viewed as a matrix, this means that 
{p^{t)*p^)/N (X, Y = N, Z, or Q) can be represented as 



F^{t) 
q 





F^'^if) F^{t) 





(5) 



where the only cross correlation is given by F^ {t) = 
{Pff{t)* p'§)/N. Thus, the dynamics associated with the 



variable Q is uncoupled from the one with N and Z, and 
is strictly connected to the single- molecule dynamics. 

Assuming the equal scattering lengths for all the con- 
stituent sites, the normalized correlator cf>^{t) is directly 
related to the cross section as measured in the coherent 
neutron scattering. The functional forms of (pg{t) and 
(/)^(t) have been chosen so that their small- wavenumber 
limits reduce to the Ist-rank reorientational correla- 
tors [23] respectively for the directions associated with 
the Y and X axes in Fig. 1. Dynamical features of 4>^{t) 
based on MD simulations have already been discussed in 
Refs. [21, 23], and some results on ^g{t) and ^^(f) in 
Ref. [23]. 

Two formulae shall be quoted here which are useful in 
analyzing simulation results for correlators. (For simplic- 
ity, we shall write down only those formulas for (pfit), 
but similar ones hold for all the correlators introduced 
above.) The correlator cf)^{t) in supercooled states ex- 
hibits the two-step relaxation: the relaxation toward the 
plateau, followed by the final relaxation from the plateau 
to zero. The latter is referred to as the a-relaxation. The 
von Schweidler law as derived from MCT including the 
next to leading order corrections [35] 



ff^-hfit/ry + h^^'Ht/rr + Oiit/rr"), (6) 



describes the departure from the plateau value - also 
referred to as the critical nonergodicity parameter - in 
the early a-relaxation region. Here and /iq denote 
the critical and correction amplitudes, respectively, b the 
von Schweidler exponent, and r the a-relaxation time. 
The exponent b is uniquely related to the so-called expo- 
nent parameter A defined in MCT [1]. Another formula, 
which is purely empirical and well adopted for fitting cor- 
relators in the a-relaxation region, is the Kohlrausch law 



(7) 



with the correlator-dependent a-relaxation time and 
the stretching exponent /3^. It is worthwhile to men- 
tion that, in the large-g limit, the Kohlrausch law (7) 
can be derived from MCT [36]. In particular, one gets 
limq^oo 0^ = b irrespective of the choice for X. Since 
the exponent b is uniquely related to A, the Kohlrausch- 
law fit in the large-g regime thus provides a means to 
estimate A based on simulation results for 4'^{t). 



III. SUMMARY OF SIMULATION RESULTS 

In this section we summarize results of MD simulations 
performed for LW OTP: more complete description of 

the simulation results can be foimd in Refs. [21 23]. We 
also present some additional results which have not been 
considered so far. At the final subsection, we discuss that 
the simulation results for LW OTP share many features 
with experimental data for real system. 
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FIG. 2: MD simulation results for the static structure factors 
at p = 2.71 molecules/nm^ and T = 190 K. (a) Solid and 
dashed lines respectively denote the static structure factor 
and the geometrical-center static structure factor S^'^. 
In this and the following figures, the arrows indicate the peak 
positions gmax (~ 14.5 nrri^"'^) in and qac (~ 9 nni^^) in 
S^*^. The inset exhibits the temperature dependence of the 
peak heights for (filled circles, left scale) and Sf^ (open 
squares, right scale), (b) Solid and dashed lines respectively 
denote the static structure factors Sq and Sq'^. Dotted lines 
refer to their self-parts, and Wq^. 



has a peak at g = qcc (~ 9 nm"-'^) which can be related 
to the inverse of the van der Waals radius rw = 0.37 nm 
for OTP molecule [37] since qqc x (2rw) ~ 7, i.e., it is 
connected to the overall size of the molecule ( cf. Fig. 1). 
Figure 2(b) exhibits the static structure factor and 

F^^{0). Also added in this 
z „„A „,,NZ Qj^g recognizes 



the cross correlation S^'^ 

figure are their self-parts, 
„„j „„z 



and 



that Sg and Wg are nearly the same. It is also seen that 
S^'^ « w^^ holds to a reasonable extent, and that the 
magnitude of the cross correlation is rather small. This 
means that, within the description based on the site-site 
static structure factors, the representation (5) for f = is 
nearly diagonal, in which essentially only is associated 
with the intermolecular correlation. 

The peak at q = qgc in S^*^ also reflects inter- 
molecular static correlations in the system. However, 
as discussed in Ref. [21] and reproduced in the inset of 
Fig. 2(a), the most pronounced temperature dependence 
in the static structure factors shows up at q = q^ax in , 
whereas the peak at g = gcc in S^'~^ is nearly tempera- 
ture independent. Thus, for LW OTP, the slowing down 
of the dynamics upon lowering T is connected with the 
intermolecular correlation manifested as the main peak 
in (the cage effect). 

In the present work, we shall primarily be interested in 
the collective dynamics arising from intermolecular cor- 
relations. As noticed in connection with Eq. (5), the cor- 
relator 4''^{t) is strictly connected to the single- molecule 
dynamics, (/)^(t) = (/'^^(i). Furthermore, « and 
shown in Fig. 2(b) imply that also 
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the correlator (f>g{t) approximately reflects the single- 
molecule dynamics only, i.e., (j)q{t) ~ (pq^s{t). Indeed, we 
confirmed both from simulation and theoretical results 
that such approximation holds to a very good extent for 
the whole time region. We shall therefore not consider 
the correlators (f)q{t) and (t)^{t) any more in the following. 



We first briefiy review the main features of the static 
structure factors, which turn out to be important in un- 
derstanding simulation results for dynamics. These static 
structure factors are also to be used as input in theoret- 
ical calculations presented in Sec. IV. 

Representative MD simulation results for the static 

structure factors in a supercooled state are presented in 
Figs. 2(a) and (b). We show in Fig. 2(a) the static struc- 
ture factor and the geometrical-center (GC) static 
structure factor 5*^^. The latter is defined by S^'^ = 
J2jj{e-'^-^^^''-^''^''>)/N with 7=;'='^ denoting the GC po- 
sition rf c = (l/3)X;Li^/ of the LW OTP molecule. 
The GC position actually coincides with the center-of- 
mass position, but we prefer to use the term "GC" be- 
cause of the reason discussed in Ref. [23] . has a main 
peak a.t q = qmax (~ 14.5 nm"^), which is compatible 
with the inverse of the LJ diameter a = 0.483 nm of a 
site in the sense that gmaxC ~ 7. On the other hand, S^'^ 



B. Density dependence of Tc and A 

Simulation results for the MCT critical temperature 
Tc and the exponent parameter A are often determined 
from the fit of the diffusion constants D according to the 
MCT asymptotic power law [1] 

D{T)^{T-T,)\ (8) 

where the exponent 7 is uniquely related to A. Such a fit 
has been performed in Ref. [22] for various densities, and 
circles in Figs. 3(a) and (b) denote the resulting Tc and A 
as a function of the density p. The critical temperature 
Tc increases with increasing p. The exponent parameter 
A seems to increase with increasing p, but the statistical 
errors in simulations do not allow to rule out the pos- 
sibility of a constant value (see below). We also notice 
that an unbiased three-parameter fit based on Eq. (8) 
suffers from correlations between the fit parameters Tc 
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Fc = 1.57 ±0.05, where the prefactor in Eq. (9) has been 

chosen to be CTpep^"' with or = 0.76 nm and er = 600 
K as in Ref. [11] for later comparison with experimental 
resuh. Thus, to a reasonable extent, the parameter F for 
LW OTP is also found to be nearly constant at the MCT 
critical point. 

Concerning the p dependence of the exponent parame- 
ter A shown in Fig. 3(b), we notice that a different value 
for A is occasionally obtained from some other analysis 
of simulation results. For example, one gets another es- 
timate for A from the Kohlrausch-law fit of some density 
correlators as described just after Eq. (7). A obtained 
in this way, based on the Kohlrausch-law fit of the cor- 
relators are also included as squares in Fig. 3(b), 
which are nearly p-independent. (We confirmed that the 
stretching exponents in the large-g regime do not depend 
on the choice of correlators.) The difference between cir- 
cles and squares in Fig. 3(b) can be considered as a sort 
of error bars in determining A based on the simulation 
results, but we already notice here that the insensitivity 
of A to density is consistent with the experimental and 
theoretical results to be described below. 
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FIG. 3: MD simulation results for (a) the MCT critical tem- 
perature Tc and (b) the exponent parameter A as a function 
of density p. Circles in (a) and (b) denote Tc and A as de- 
termined from the fit of the diffusion constants T) according 
to Eq. (8). In (b), squares refer to A based on the expo- 
nent 6, which is determined from the largc-g behavior of the 
Kohlrausch stretching exponents \\va.q^oo Pq=b for the den- 
sity correlators 4'^{t). In both (a) and (b), lines are guide to 
the eyes. 

and 7 [14]. It also suffers from the identification of a 
fitting T range bounded below from Tc and above from 
the T at which correlation functions start to obey the 
time-temperature superposition. We note in passing that 
deviation from the MCT behavior for T < Tc has been 
discussed for the LW-OTP model in Ref. [22]. 

Often, the p dependence of Tc can be condensed into 
an effective coupling parameter 

F oc pT-i/^ . (9) 

This is the only relevant parameter for specifying the 
thermodynamic state of soft-sphere systems whose re- 
pulsive interaction is proportional to r~^^ [38]. Indeed, 
in the case of binary mixture of soft spheres, it was found 
from computer simulations that the ideal glass transition 
occurs at a constant value, F = Fc [39, 40]. Although 
Eq. (9) is valid only for the r~^'^ soft-spliere system, a 
computer-simulation study for a model of polymer at dif- 
ferent pressures [41] indicate that, also for LJ systems, 
the MCT critical point might be described well in terms 
of the effective coupling parameter. For LW OTP, we 
find from the p dependence of Tc shown in Fig. 3(a) that 



C. Critical nonergodicity parameters and critical 
amplitudes 

Figures 4(a) and (b) show the q dependence of the crit- 
ical nonergodicity parameters /^'^ and the critical ampli- 
tudes of the correlators (t) based on the fit accord- 
ing to Eq. (6) for three representative densities, p = 2.71, 
2.83, and 2.97 molecules/nm^. The fit of (t)^{t) accord- 
ing to Eq. (6) has been performed by constraining the 
exponent b to the value specified by A, the latter being 
taken from the squares in Fig. 3(b), and by regarding f^'^, 

/t'^ and hq^'^^ /t'^^ as fitting parameters. Thus, and 

h^^"^"^ can be determined only up to (j-independent mul- 
tiplicative factors, and this is why shown in Fig. 4(b) 
are given in arbitrary units. (Hence, one cannot directly 
compare the amplitude but only the wavevector depen- 
dence in for different densities.) 

For q > ^max (~ 14.5 nm~^ for all the densities consid- 
ered) , for the three densities are very similar to each 
other. For q < gmax, on the other hand, fg'^ is higher for 
lower density, and this holds in particular around the in- 
termediate wavenumber qqc (~ 9 nm"-'^): whereas f^" at 
q ~ qcc exhibits only a shoulder for the highest density 
p = 2.97 molecules/nm^, a well-developed peak is dis- 
cernible for the lowest density p = 2.71 molccules/nm'^. 

The results for shown in Fig. 4(b) are somewhat 
noisy compared to /^^. Concerning the q dependence, 
we notice that there are two minima in located at 
Q ~ Qgc and (/max- As will be mentioned in Sec. Ill D, this 
is consistent with the existence of the peaks in f^'^ around 
these two wavenumbers. Furthermore, the minimum in 
/ig at g ~ qcG is more pronounced for lower density in 



6 



1 

0.8 
0.6 

0.4 
0.2 
1 

0.8 
0.6 

0.4 
0.2 




' 1 


1 

^ Nc 

q 


1 ' 

(a) - 


I \^ 

IgC ^mai 

" 11 




- 

1 r" — ' 


1 1 


— 1 

q 


(b) - 










\ 

"\ r\ 




1 




1 



20 40 60 
q [nm"'] 




FIG. 4: MD simulation results for (a) the critical nonergod- 
icity parameters f^'^ and (b) the critical amplitudes ftg of 
the correlators 0^ {€) as determined from the fit according to 
Eq. (6) for three densities p = 2.71 (solid lines), 2.83 (dashed 
lines), and 2.97 molecules/nm^ (dotted lines), are in ar- 
bitrary units. 



the sense that h^^^^/h^^^ is larger for lower p, and this 
is also consistent with the found density dependence of 
f^'' around this wavenumber. 



FIG. 5: MOT results for the wavenumber dependence of var- 
ious quantities for the hard-sphere system as a function of 

qd with d denoting the hard-sphere diameter, (a) Solid and 
dashed lines respectively denote the critical nonergodicity pa- 
rameters fq and the critical amplitudes hq for normalized 
density correlators 4>q{t). Dotted line refers to the Percus- 
Yevick static structure factor Sq, multiplied by a factor of 
0.08 for ease of comparison, at the critical point, (b) Solid 
and dashed linos respectively denote the a-rclaxation times 
Tq and the stretching exponents Pq obtained from the fit ac- 
cording to Eq. (7) of the a- master curve for 4>q{t). Tq are in 
arbitrary units. 



D. Unusual feature in the wavenumber dependence 

All parameters describing density correlations exhibit 
a characteristic wavenumber dependence reflecting that 
of the underlying static structure factor. As found in 
Ref. [21] and also discussed in more detail in Ref. [23], 
however, the dynamics in LW OTP exhibits an unusual 
wavenumber dependence of the critical nonergodicity pa- 
rameters and the a-relaxation times of the correlators 
(j)^ (t) at intermediate wavcnumbcrs q qqc ■ Before cm- 
barking on the unusual feature, let us first consider a 
"usual" case as a reference. 

We show in Figs. 5(a) and (b) the critical nonergod- 
icity parameters /°, the critical amplitudes hq, the a- 
relaxation times Tq, and the Kohlrausch stretching expo- 
nents /?q of normalized density correlators <j)q{t) for the 
hard-sphere system as determined by solving the MCT 
equations for simple systems [1, 42]. Tq and (3q have been 
obtained from the fit according to Eq. (7) of the a-master 
curve for (pq{t). The static structure factor used is eval- 
uated within the Percus-Yevick approximation [38], and 
the one at the critical point is included in Fig. 5(a). A 



strong correlation in the q dependence of these quanti- 
ties can clearly be observed for the whole wavenumber 
regime: 1/hq, Tq, and (Jq oscillate in phase with Sq, 
and this holds in particular around the first sharp diffrac- 
tion peak, qd ^ 7 with d denoting the hard-sphere diam- 
eter. A theoretical explanation of such correlation has 
already been documented [29, 35], and shall not be re- 
peated here. Such "usual" results have been observed in 
simulation studies for LJ binary mixture [43], silica [44], 
and water [45], and also in experimental results reviewed 
in Ref. [46]. 

The result found in the; sinnilation study for LW 
OTP is unusual in that such a correlation is violated 
at intermediate wavenumbers q « gcc- This is sum- 
marized in Figs. 6(a) and (b) for the density p = 
2.71 molecules/ nm"^. A related figure for p = 2.83 
molecules/nm^ can be found in Ref. [23]. One recognizes 
from Figs. 6(a) and (b) that f^", and of the 

correlator ^^(t) exhibit an additional peak at g « ^gc 
which does not exist in the corresponding static struc- 
ture factor . A signature of such a peak in is also 
discernible, but a definitive conclusion cannot be drawn 
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FIG. 6: MD simulation results for the wavenumber de- 
pendence of various quantities for LW OTP at p = 2.71 
molecules/nm^. (a) Solid and dashed lines respectively de- 
note the critical nonergodicity parameters fg'^ and the criti- 
cal amplitudes for the correlators (J)^{t). arc in arbi- 
trary units. Dotted line refers to the static structure factor 
, multiplied by a factor of 0.08 for ease of comparison, at 
T = 190 K. (c/. Tc w 172 K for this density, see Fig. 3.) (b) 
Solid and dashed lines respectively denote the a-relaxation 
times Tq and the stretching exponents /3q obtained from the 
fit according to Eq. (7) of the correlators 0^(t) at T = 190 K. 

since the results for (3^ are rather noisy. As already no- 
ticed in connection with Figs. 4(a) and (b), this unusual 
feature is found to be more pronounced for lower density. 
It was suggested in Refs. [21, 23] that the unusual fea- 
ture is caused by the coupling of the rotational motion 
to the GC motion. This follows from the fact that the 
GC static structure factor S^^ has a peak at g = qoc as 
shown in Fig. 2(a). 

To understand what is really the unusual feature, 
especially concerning f^" shown in Fig. 6(a), we go 

back to the original site-site density correlators F^^it). 
Let us recall that (l)^{t) = F^{t)/S^ with F^(t) = 
Ea,6=i^^f (i)/3 and S^^ = Ea,i,=i 5f /3. Critical non- 
ergodicity parameters F^^'^ of the correlators Fg^{t), ob- 
tained similarly from the fit according to Eq. (6), and 
the site-site static structure factors S^^ are shown in 
Figs. 7(a) and (b), respectively. (Only the independent 
components in F^^'"^ and S^^ arc shown. The indepen- 
dent components, e.g., in the total nine 5°'', are Sg^, 
, Sf, and S'f due to the symmetry Sf = 5^" as 
well as the C2v symmetry of the LW OTP molecule.) It 
is seen from the comparison of Figs. 7(a) and (b) that the 



FIG. 7: MD simulation results for the site-site critical non- 
ergodicity parameters F"*"^ and the site-site static structure 

factors Sq'' for the density p = 2.71 molcculcs/nm''. (a) Fq'^" 
(solid line), F^'^" (dotted line), F^'^" (dashed line), and f}^" 
(long-dashed hne). (b) 5" (solid line), 5,^ (dotted hue), 5,^ 
(dashed hne), and S^^ (long-dashed hne) at T = 190 K. 



results for LW OTP are "usual" in the sense that the q 
dependence of F^^ " correlates well with that of S^^ in the 
whole wavenumber regime including q « qcc- We also 
notice that these quantities take positive and negative 
values at g « qcc- The unusual feature concerning f^" 
shows up only after summing up the site-site F^^'^ and 
Sg'' to obtain F^'^ and 5*^. One observes that positive 
and negative components in S^^ almost cancel out after 
taking the summation, which results in small and flat 5^ 
at q K qQC as shown in Fig. 6(a). On the other hand, 
such a cancellation does not occur in the summation of 
the components i^"*"^, and this causes the unusual peak 
in and hence in /^^ = F^^/S^, at 9 « gcc- Thus, 
the mentioned unusual feature reflects purely dynamical 
effects, since it can be observed only in such dynamical 
quantities as f^", l/h^ and . 



E. Comparison with experimental results 

Sine the LW OTP model is a very simplified one, one 
might think that a straightforward comparison of the MD 
simulation results with experimental data is not feasible. 
However, we discuss in the following that the simulation 
results for LW OTP share many features with experi- 
mental ones. Since we are interested in the q dependence 
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of collective dynamical quantities, especially in the vicin- 
ity of qcc and gmax, we shall be mostly concerned with 
coherent neutron-scattering results. A review of neutron- 
scattering studies of OTP is presented in Ref. [5] , and the 
experimental results to be presented below can be found 
in this article and references cited therein. 

Let us first consider the static structure factor Sg^^ 
from coherent neutron scattering for fully dcuterated 
OTP, which is given as a weighted sum of atomic cor- 
relations, Sg'^^ oc J2a,0 babfsSg^. Here a and /3 refer 
to deuteron and carbon atoms, ba the scattering length, 
and Sg^ the partial structure factors. (Greek characters 
a and /3 are used here to distinguish them from Roman 
characters a and b which have been adopted to label sites 
in the LW OTP molecule. Also, quantities as determined 
from experiments shall be distinguished with the super- 
script "exp".) It has been observed that, in contrast 
to atomic systems, the main peak of S^^^ is split into 
two maxima at about g = 14 and 19 nm~^ [4]. Also, 
a shallow shoulder is discernible in Sg^^ around q = 9 
nm^^. It was conjectured that the peak at 17 » 19 nm~^ 
is built up mainly by intramolecular correlations within 
phenyl rings, while the maximum at g ~ 14 nm~^ is as- 
sociated with intermolecular correlations between phenyl 
rings. The shoulder at g « 9 nm^^ was interpreted as 
being due to correlations between molecular centers of 
mass since its position is compatible with the inverse of 
the van der Waals radius rw = 0.37 nm [37]. 

This picture for S^'^^ is consistent with the simulation 

results for and S^''^ shown in Fig. 2(a). Since each 
site of the LW OTP molcciilc represents an entire phenyl 
ring, the peak at gmax ~ 14.5 nm"-'^ in corresponds 
to the first main peak in Sg^^. No second peak aroimd 
g w 19 nm^^ can be observed in since the internal 
structure within a phenyl ring is completely discarded in 
the LW OTP model. The peak at gcc ~ 9 nm"! in S^'^ 
is related to the shoulder in Sg^^, although this is not re- 
flected in . As discussed in connection with Fig. 7(b), 
the disappearance of any shoulder or peak at g « gcc in 
is due to the cancellation of the constituent site-site 
correlation functions Sg^, which do exhibit (positive and 
negative) peaks around this wavenumber. Such an almost 
perfect cancellation might not occur in Sg^^, due to more 
complicated nature of the constituent atomic correlations 
Sg^ in real system. Indeed, this conjecture is supported 
by MD-simulation studies for more realistic OTP mod- 
els [47, 48], where a shoulder at g w 9 nm^^ is discernible 
in the static structure factor which corresponds to Sg^'^. 

It has been discussed for a liquid of linear molecules 
that orientational correlations can lead to a prepeak at 
low g [49]. By calculating the corresponding static orien- 
tational correlation functions, we confirmed that there is 

no prominent peak at gcc ~ 9 nm"-'^ in these functions 
for LW OTP. Therefore, we do not think the shoulder at 
g « 9 nm~^ as observed in Sg^^ reflects the orientational 
correlations of the kind discussed in Ref. [49] , and this is 
consistent with the picture that the shoulder stems from 



the correlations between the molecular centers of mass. 

We next compare the simulation and experimental re- 
sults for the p dependence of Tc and A. It has been shown 
in Ref. [11] from an analysis of incoherent density correla- 
tors at various pressures that the p dependence of Tc can 
be combined to the effective coupling parameter given in 
Eq. (9). This observation is in agreement with the simu- 
lation result for LW OTP. Furthermore, using the same 
prcfactor as described just after Eq. (9), the experimen- 
tal value r°^P f« 1.498 ± 0.004 characterizing the MCT 
critical point is rather close to Tc ~ 1.57 ± 0.05 found 
for LW OTP. A smaller error bar in F^p might be due 
to a narrower density range investigated in the experi- 
ment. Concerning A, its insensitivity to density has been 
demonstrated with the experimental value A*"^? « 0.77. 
This is in accord with the simulation result shown as 
squares in Fig. 3(b), including the value for A. 

Coherent as well as incoherent neutron-scattering re- 
sults for density correlators in supercooled states exhibit 
two-step relaxation in agreement with the prediction of 
MCT. The wavenumber dependence of various quantities 
characterizing such glassy dynamics has been determined 
from fits of those density correlators according to MCT 
asymptotic formulas or to the Kohlrausch law, like we did 
for simulation results. It was observed that the critical 
nonergodicity parameters fg^'^'^ of coherent density cor- 
relators oscillate in phase with S"^"^, with the two peaks 
around gcc and gmax [4]. The critical amplitudes /i^^p 
were found to oscillate in antiphase with S'^'^p, with the 
existence of two minima at g ~ gcc and gmax- These 
trends are in agreement with the simulation results for 
0^(t) shown in Fig. 6(a). 

Let us now consider the g dependence of the relaxation 
times and the stretching exponents of coherent density 
correlators in the a-relaxation regime. Some reservation 
is necessary in the experimental results for small g due 
to the presence of incoherent background and to contri- 
butions from multiple scattering. Therefore, the exper- 
imental a-relaxation times r|^P show a tendency to in- 
crease for decreasing g. Nevertheless, on top of such back- 
ground, it was observed that t^'^p exhibit two plateaus 
around gcc and gmax [12]. This is a signature of the 
characteristic g dependence as found in for LW OTP 
shown in Fig. 6(b). Concerning the stretching exponents 
(3'^^, a systematic variation in phase with S'^'^p was ob- 
served, especially in the vicinity of gmax [12], although a 
definitive conclusion cannot be drawn for g « goc due to 
the noise in the experimental results. Thus, the overall 
g dependence of t^'^p and jS"^"^ is in accord with the one 
for LW OTP shown in Fig. 6(b), and even a signature of 
the unusual peak at g « gcc as discussed for LW OTP 
can be observed in Tg^'^. 

A natural question arises as to whether the peaks at 
g ~ Qgc found in the experimental results for /^p'^, 
1 / /iq'^P and Tg'^P can really be considered as unusual. This 
is because, in contrast to for LW OTP, the experi- 
mental static structure factor Sg^^ exhibits a shoulder, 
though tiny, in this wavenumber regime. However, it is 
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rather surprising that such a tiny shoulder in S'^p is re- 
lated to a pronounced wavenumber variation of those dy- 
namical quantities, and it seems worthwhile to pay spe- 
cial attention to dynamical features around qcc- There- 
fore, we consider that a further investigation for the un- 
usual features in LW OTP is valuable and might also be 
relevant in understanding the experimental results. 

IV. THEORETICAL RESULTS 

In this section, results for dynamical quantities cal- 
culated from MCT based on the site representation are 
presented, and are compared with those from the MD 
simulations. In particular, we examine whether the un- 
usual feature at intermediate wavenumbers discussed in 
Sec. Ill D can be accounted for by the theory. Since the 
unusual feature is found to be more pronounced for lower 
density, most of the theoretical calculations shall be done 
for the lowest density p = 2.71 molecules/nm^ studied in 
the MD simulations. 

A. MCT equations based on the site representation 

Within the site representation for molecules, the dy- 
namics of the system is most naturally characterized 
by the site-site density correlators F^^{t) defined for 
a, 6 = 1, • • • , n, where n denotes the number of sites in 
a molecule. The MCT equations for F^^{t) consist of an 
exact Zwanzig-Mori equation and an approximate expres- 
sion for the relaxation kernel, which are given in Ref . [33] . 
Regarding F^^{t) as elements of an n x n matrix Fq(t), 
the Zwanzig-Mori equation is given by 

d'l¥q{t) + ¥^{t) + f dt' mgit - t') dt'Fg{t') = 0, 

Jo 

(10a) 

where denotes the characteristic frequency matrix 

nl = q^J,S-\ (10b) 

with S^"^ representing the inverse matrix of Sq. ./^^ de- 
note the site-site static (longitudinal) current correlation 
functions, whose explicit expressions for molecules pos- 
sessing the C2v symmetry (like LW OTP and water) in 
terms of molecule's inertia parameters can be found in 
Ref. [50]. The MCT expression for the relaxation kernel 
reads 

Ulgit) = SgJ^g[Fit)], (lla) 

where the mode-coupling functional J^q is given by the 
equilibrium quantities: 

<A>'(g;fc,p) = (^{«"*- [^"^'fccf + 5«Vcf' ]} 

X {g . [5'>^'kc{^ + (5''V4^']}/g^ (He) 



with p = q — k. (We have used here a slightly differ- 
ent convention for writing the mode-coupling functional 
from the one adopted in Ref. [33] to simplify some equa- 
tions which follow.) Here, the direct correlation function 
is defined via the site-site Ornstein-Zernike equation [38] , 
pcf = [w-^]"'' - [S-i]"''. Equations (10) and (11) con- 
stitute a set of closed equations of motion for determin- 
ing Fq{t), provided the static structure factors S^'' are 
known. In the present work, S^^ determined from MD 
simulations shall be used. 

The matrix of long-time limits (or the nonergodicity 
parameters), = Fg(t ^ oo), obeys the implicit equa- 
tion defined by the mode-coupling functional J-g, 

Fq[Sq-Fq]-'=SqJ^q[F]. (12) 

This equation can be derived from Eqs. (10a) and (Ha) 
by taking the t ^ oo limit. From an iterative proce- 
dure F^^+^^[S, - F^^'+^']-^ = S9:F,[F(j)] starting with 
Fg°^ = Sq, one obtains a solution of Eq. (12) as Fg = 

limj^oo Fq"'' [33]. One gets trivial solutions F^ = for 
T > Tc, whereas nontrivial solutions Fg ;^ can be ob- 
tained for T < Tc- Here, and in the following, we mean 
by F, ;^ (or F^'' >- 0) that the matrix Fg is positive 
definite. Thus, one obtains as the highest temperature 
at which there holds Fg 0, and the solution at this crit- 
ical point provides the critical nonergodicity parameters 

F". 
q- 

The convergence of the iterative procedure for Eq. (12) 

is ruled by the spectral radius of a stability matrix, which 
can be defined from the mode-coupling functional J-g as 
in the case of simple systems [1] and is given by 

p \,fl,X' ,IJI,' 

X [Sfc - Ffe]^«' [Sfe - Fkf'' F^''^' . (13) 

In deriving this expression from Eq. (lib), the wavevec- 

tor integrals arc converted into discrete sums by intro- 
ducing some upper cutoff and using a grid of M equally 
spaced values for the wavenumbers. Thus, the wavenum- 
ber (g, k and p) can now be considered as a label for 
an array of M values. Correspondingly, the coefficients 

Vx^X'A^'^ ^'^) (ll'^) expressed as V;a%;t.' ; 

in Eq. (13). The details of the transformation of the 
mode-coupling functional to a polynomial in the dis- 
cretized variables can be found in Ref. [35]. 

For the notational simplicity, we introduce new indices 
i = {q,a,b) and j = {k,a',b') using the so-called dictio- 
nary order, which run from 1 to Mn^. Then, the stability 
matrix given in Eq. (13) is simply denoted as Cjj. Unlike 
the case of simple systems for which the stability matrix 
is given by positive matrix [1], each element of Cij can 
take positive and negative values. However, can be 
considered as a generalized positive matrix in the sense 
that it transforms a positive definite matrix into another 
one: if a;i ;^ (meaning x^'' >- 0), then J2j ^ij^j ^ 
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and XiCij >~ 0. Such matrix Cij has a nondcgcncratc 
maximum oigcnvahic E < 1, and the corresponding right 
(ej) and left (ci) eigenvectors can be chosen as ;^ and 
ej !^ [51]. The MCT critical point is characterized by 
= 1. Let Ci and it specifically denote the right and 
left eigenvectors, respectively, of the stability matrix Cfj 
at the critical point: ^fj^j = J2i ^i^fj — ^i- "^^^ 
eigenvectors are fixed uniquely by requiring CiCi = 1 
and ti [e(S — F°)e]i = 1. These eigenvectors can be 
used evaluate the critical amplitudes 



q — I'-'q - Fg] l^q " ^g] > 



H„ = [S, 

and the exponent parameter A 



A 



EE 

q a,b 



q 



ah 



HI 



(14) 



(15) 



The above formulation for molecular systems is essen- 
tially the same as the one for simple systems, the only 
difference being the appearance of matrices in place of 
scalar quantities. It is then obvious that all the uni- 
versal results concerning the MCT-liquid-glass-transition 
dynamics, as originally derived for simple systems [1], are 
valid also for MCT for molecular systems with such re- 
placement of scalar quantities with matrices. 

Since the LW OTP moleeide consists of three sites, 
a natural choice would be n = 3, and the MCT equa- 
tions (10) and (11) for this model are given by 3 x 3 
matrix equations. On the other hand, the analysis of the 
simulation results presented in Sec. Ill D concerning the 
unusual feature at intermediate q « gcc suggests the im- 
portance of taking into account the spatial correlation of 
GC through the static structure factor S^'^ . This implies 
that, in accounting for the unusual feature at g w qGCi 
it might be necessary to include GC as the additional 
4th site in the theoretical calculations. This leads to an- 
other formulation for LW OTP based on 4 x 4-matrix 
site-site density correlators i^^''(i), for which the MCT 
equations (10) and (11) are given by 4 x 4 matrix equa- 
tions. The required static inputs S'^^ for a, 6 = 1, • • • , 4 
in the new formulation can also be determined from MD 
simulations. By comparing theoretical results with and 
without including GC, one can judge whether the un- 
usual feature at the intermediate wavenumbers q ~ qcc 
discussed in Sec. Ill D stems from the coupling to the GC 
dynamics or not. 

One might think that the inclusion of GC as an ad- 
ditional site is an ad hoc procedure, which is motivated 
in view of the simulation results. This kind of problems 
can occur in the site-density formulation since the site- 
density fluctuations pi defined for finite number of sites 
do not provide a complete set of variables describing the 
dynamics of molecules. This is in contrast to the tensor- 
density fluctuations as adopted in Refs. [17 19, 24 29], 
which do provide a complete set of variables for molecules 
and also naturally incorporate the GC correlation. How- 
ever, as mentioned in Sec. Ill A, the most important den- 
sity fluctuations relevant for the structural slowing down 
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FIG. 8: Same as in Fig. 3, but here theoretical results are 
included as well. In both (a) and (b), MCT results with and 
without including GC are denoted as filled diamonds (con- 
nected with thick solid line) and filled triangles (connected 
with thick dashed line), respectively. 



manifest themselves as the main peak at q = (/max in 
(the cage effect). Such correlations arc incorporated in 
the site-density formulation even with the natural choice 
for the number of sites (n = 3 for LW OTP). Indeed, it 
will be shown below that the MCT equations with n = 3 
can describe the basic feature of the simulation results, 
and this formulation is already useful. Our attitude here 
to include GC as an additional site (resulting in n = 4) 
is only for the investigation of the dynamical features at 
q ~ qcc- In fact, we will see that the inclusion of GC 
does not alter significantly the results for the wavenum- 
bers other than q w qcc- 



B. Density dependence of Tc and A 

Figures 8(a) and (b) respectively show the theoreti- 
cal results for and A as a function of the density p, 
along with the corresponding simulation results to be de- 
noted as and A''^'^ from here on in this section. It 
is seen from Fig. 8(a) that both of the theoretical Tc, 
with and without including GC, exhibit qualitatively the 
same density dependence as that of T^°. Concerning 
A shown in Fig. 8(b), the theoretical results with and 
without including GC are practically the same, and are 
nearly independent of p. The latter feature is in accord 
with the simulation result shown as squares in Fig. 8(b). 
The agreement between the theoretical and simulation 
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results for the value of A is reasonable in view of the er- 
ror bars in estimating A'^^ as discussed in connection 
with Fig. 3(b). 

The theoretical Tc with including GC is in better agree- 
ment with T^^, but is still located considerably below 
at all the densities investigated. One might think 
that the discrepancy between the theoretical and 
is much larger than the one known for the hard-sphere 
system (HSS): MCT for HSS yields the critical packing 
fraction ipc which differs only about 7% from the experi- 
mental value (^^'^P for hard-sphere colloids [14], for which 
HSS is known to serve as a good model. (The previous 
estimate of the difference between the theoretical ipc and 
the experimental (/p^'^p was about 15% [1], but a recent 
analysis performed in Ref. [14] indicates that half of this 
error is due to the use of the Percus-Ycvick static struc- 
ture factor as input instead of simulated one.) However, 
we argue in the following that the discrepancy between 
the theoretical and T^^ for LW OTP is of comparable 
size to the one for HSS. 

To this end, we notice that the p dependence of 
could be condensed into a nearly constant effective cou- 
pling parameter Vf^ = 1.57 ± 0.05 (c/. Sec. IHB). We 
found that, to a reasonable extent, this holds also for the 
theoretical results: from the p dependence of the theo- 
retical Tc, one gets T,, = 1.77 ± 0.07 and 1.9 ± 0.1 with 
and without including GC, respectively. Thus, in terms 
of the critical effective coupling parameter Fg, the dis- 
crepancies between the theoretical and simulation results 
are only within 13% and 20% with and without includ- 
ing GC, and arc of comparable size to the discrepancy 
found for HSS. A similar analysis has been done for bi- 
nary mixture of LJ particles (BMLJ) [15], for which it 
was discussed that the difference between Tg = 0.922 
from MCT and T^^ = 0.435 from the MD simulation 
is comparable (w 20 %) to the one for HSS when quan- 
tified in terms of the effective coupling parameter. We 
thus conclude that our theoretical estimate of for LW 
OTP is within the comparable error bar as the one for 
HSS and BML.J. 

However, we notice that our theoretical results exhibit 
an unconventional feature in that is underestimated 
compared to T^^ . This is in contrast to all previous 
MCT calculations, where MCT is found to overestimate 
Tc (or underestimate when the packing fraction is con- 
cerned): e.g., ifc = 0.546 < VJ^^P fa 0.58 for HSS [14], and 
Tc = 0.922 > T^^ = 0.435 for BMLJ [15]. In particu- 
lar, this is also in contrast to MCT results for molecules 
based on the tensor-density formulation [17-19]. Con- 
cerning this problem, we notice that including the triple 
direct correlation function C3 was found to move up T,, for 
LW OTP more than by a factor of 2 [21]. Although the 
finding in Ref. [21] is based on the simplified scalar MCT 
equations dealing with the correlators 0^ (t) only, we ex- 
pect that the inclusion of C3 would significantly move up 
Tc shown in Fig. 8(a) and would lead to Tc > T^^ in 
agreement with the previous MCT studies. 

Such theoretical calculations with including C3, how- 
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FIG. 9: Comparison of MD simulation and theoretical results 
for (a) the critical nonergodicity parameters f^'^ and (b) the 
critical amplitudes of the correlators <t>^ (t) for the density 
p = 2.71 molecules/nm^. In both (a) and (b), circles denote 
the simulation results as determined from the fit according 
to Eq. (6), while solid and dashed lines refer to MCT results 
with and without including GC, respectively. The simulation 
results for hq are in arbitrary units. 



ever, shall not be performed in the present work. This is 
basically because an accurate evaluation of C3 from MD 
simulations is qiutc a demanding task. Furthermore, it 
was also found in Ref. [21] for LW OTP that including C3 
does not significantly alter theoretical results other than 
Tc, and we consider that MCT without C3 captures the es- 
sential physics in the dynamics of supercooled LW OTP. 
Indeed, as we will see below, our MCT results without 
the use of C3 are in semiquantitative agreement with the 
simulation results. 



C. Critical nonergodicity parameters and critical 
amplitudes 

Figures 9(a) and (b) compare the theoretical and MD 
simulation results for the critical nonergodicity parame- 
ters /^'^ and the critical amplitudes of the correlators 
4>^{t), respectively. In both of these figures, solid and 
dashed lines respectively refer to the MCT results with 
and without including GC, while circles denote the sim- 
ulation results. 

Concerning f^" shown in Fig. 9(a), it is seen that, even 
without including GC, the q dependence of the simula- 
tion result is well reproduced by the theory, especially 
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for wavenumbcrs q ^ (/max- However, the theory without 
GC fails to reproduce the peak in f^"^ at q ^ ^gc found 
in the simulation result. This peak is reproduced by the 
theory which includes GC, although its magnitude is un- 
derestimated. Furthermore, the overall agreement with 
the simulation result becomes better also for q > g'max by 
including GC. 

Also for exhibited in Fig. 9(b), the overall q depen- 
dence of the simulation result is well reproduced even by 
the theory without GC. On the other hand, the mini- 
mum in a,t q Ki groc found in the simulation result is 
accounted for only by the theory which includes GC. The 
found improvement in the theoretical results for /^"^ and 
a,t q Ki qcc, which is achieved by including GC, sup- 
ports the idea that the unusual wavenumber dependence 
discussed in Sec. HID is basically due to the coupling to 
the GC dynamics. A further support shall be discussed 
in Sec. IV D from the analysis of the a-relaxation times. 



D. Dynamics in the a-relaxation region 

As discussed in Rcf. [23] for LW OTP based on the 
MD simulation, the most faithful tests of MCT concern- 
ing the dynamics should be performed in the a-relaxation 
part starting from the plateau regime. This is because 
the approach toward the plateau, for which MCT pre- 
dicts an asymptotic power-law ~ (0 < o < 0.5) [1], 
was found to be almost completely masked by the micro- 
scopic dynamics. In view of this, tests of the theoretical 
results for density correlators shall be done in terms of 
MCT a-master curves. Of particular relevance here is 
the MCT second scaling law - also referred to as the su- 
perposition principle - which states that correlators in 
the a-relaxation region for different temperatures can be 
superposed on top of each other simply by rescaling the 
time scale: 

<l>^{t) = ^fit/rf) . (16) 

Here denotes the a-master function. 

In the strict test of the MCT second scaling law, on 
the other hand, one cannot freely choose the scale t^, 
which can depend on the choice of the variable X as well 
as on the wavenumber q. According to MCT, there exists 
a single time scale, say t, characterizing the a relaxation 
of all the correlators [1]. Thus, instead of Eq. (16), one 
actually has for the MCT second scaling law 

cj^^it)^^fit/T). (17) 

The MCT a-master function <j)^(i) can be evaluated from 
the MCT equations at T = Tc up to an overall time scale, 
with the initial behavior given by the von Schweidler law, 
Eq. (6) [1]. 

The mentioned second scaling law of MCT implies 

that the test of the MCT a-mastcr functions against 
simulation results for density correlators should be done 
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FIG. 10: Comparison of MD simulation and theoretical re- 
sults for the correlators 4>q (t) for the density p = 2.71 
molecules/nm'^ at three wavenumbers q = 13.1, 20.0, and 
30.0 nm~^. Circles denote the simulation results at T = 190 
K. Solid lines refer to the a-master curves ^^{t/r) from MCT 
including GC, where r has been chosen so that both the simu- 
lation and theoretical curves at g = 13.1 nm~^ yield the same 
a-relaxation time when fitted with Eq. (7). 



by adjusting the single time scale r only. Such test is 

performed in Fig. 10 for the correlators <^^(i) at three 
wavenumbers q = 13.1, 20.0, and 30.0 nm~^, which are 
close to the first peak, the first minimum, and the second 
minimum in , respectively [cf. Fig. 2). In Fig. 10, cir- 
cles refer to the simulation results for (t>^{t) at p = 2.71 
molecules/nm^ and T = 190 K. Solid lines denote the 
a-master curves cj)^{t/T) from MCT which includes GC, 
where r has been chosen so that both the theoretical and 
simulation curves at g = 13.1 nm~^ yield the same a- 
relaxation time when fitted with Eq. (7). It is seen 
from Fig. 10 that the MCT a-master curves describe 
well the time dependence of the simulation results for 
(j)^ (t) in the a-relaxation region, including the relative 
a-relaxation times for different wavenuinl)ers. (Notice 
that, from the construction of r, the real test of the rel- 
ative a-relaxation times is performed only for q = 20.0 
and 30.0 nm~^ in Fig. 10.) In particular, the relaxation 
stretching, which is pronounced for q — 20.0 and 30.0 
nm~^, is well reproduced by the theory. Indeed, from 
the Kohlrausch-law fit according to Eq. (7), we found, 
e.g., for q = 20.0 nm-\ (3^ ^ 0.62 from the MCT a- 
master function and 0.60 from the simulation curve. A 
similar quantitative test of MCT against simulation re- 
sults for density corrc^lators, which uses the single time 
scale as an adjustable parameter, has been performed in 
Ref. [52] for binary mixture of LJ particles. 

Theoretical and simulation results for the a-relaxation 
times and the stretching exponents of the corre- 
lators (t) for the whole wavenumbers are compared in 
Figs. 11(a) and (b), respectively, which are obtained from 
the fits according to Eq. (7). In these figures, solid and 
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FIG. 11: Comparison of MD simulation and theoretical re- 
sults for (a) the a-relaxation times and (b) the stretching 
exponents of the correlators (p^ (t) for the density p = 2.71 
molcculcs/nm'*. In both (a) and (b), circles denote the simu- 
lations results at r = 190 K, while solid and dashed linos refer 
to results from the MCT a-master curves with and without 
including GC, respectively. The overall time scale of the the- 



oretical results for has been chosen so as to reproduce the 
same at g = 13.1 nm"'^ as that from the MD simulation. 



dashed lines respectively denote the MCT results with 
and without including GC, whereas circles refer to the 
simulation results. Again, the overall scale of the theo- 
retical results for the a-relaxation times has been chosen 
so as to reproduce the same at g = 13.1 nm^^ as 
that from the MD simulation. It is seen that, even with- 
out including GC, the q dependence of the simulation 
results for and fi^ is well reproduced by the theory at 
the semiquantitative level, especially for the wavenum- 



bers q 



> 



For 



q < qn 



on the other hand, it 



is seen from Fig. 11(a) that the unusual peak in the a- 
relaxation times at intermediate q w Qgc, a-s observed 
in the simulation result, is reproduced only by the the- 
ory which includes GC, although its magnitude is still 
underestimated. This again supports the idea that the 
unusual peak is basically due to the coupling to the GC 
dynamics. 



V. SUMMARY AND CONCLUDING REMARKS 

In this paper, we reported MD simulation results per- 
formed for a model of molecular liquid OTP developed 



by Lewis and Wahnstrom, paying special attention to the 
wavenumber dependence of the structural a relaxation 
of the collective dynamics, and showed that the simu- 
lation results for the model share many features with 
experimental data for real system (Sec. 111). We then 
demonstrated that theoretical results based on MCT cap- 
tures the simulation results at the semiquantitative level 
(Sec. IV): it is found that MCT yields a fair estimate of 
the critical temperature Tc and the exponent parameter 
A including their density dependence, and predicts the 
wavenumber dependence of dynamical quantities rather 
well, in particular near the first sharp diffraction peak 
9max of the static structure factor . Through these 
investigations, we established the relevance of our theo- 
retical results and their interpretation in understanding 
experimental data for real system. 

As described in Sec. Ill A, major intermolecular cor- 
relations manifest themselves as the peaks at q = ^max 
in and at g = qcc in S^'"- On the other hand, the 
most pronounced temperature dependence in the static 
structure factors shows up at g = gmax in {cf. the 
inset of Fig. 2(a)). This indicates that the structural 
slowing down and anomalous glassy features in the dy- 
namics upon lowering T are primarily caused by the in- 
termolecular correlation manifested as the main peak in 
(the cage effect). This is supported by the observa- 
tion in Sec. IV that the theoretical results, which do not 
take into account S^'^ , already capture the basic features 
of the simulation results. 

On the other hand, though of subordinate nature in 
the above sense, the simulation results for LW OTP ex- 
hibit interesting and unusual properties at intermediate 
wavenumbers q « qcc, which reflect purely dynamical ef- 
fects. As discussed in Sec. HIE, similar features can also 
be observed in experimental data for real system. We ar- 
gued that such unusual features for q « gcc are basically 
due to the coupling to the GC dynamics. This is because, 
compared to the simulation results, the theoretical results 
for q « qcc were found to be improved by including the 
spatial correlation of GC through S^''^ . However, there 
still remain quantitative discrepancies between the theo- 
retical and simulation results for q w qcc compared with 
the agreement found for the other wavenumber regime. 
This implies that the present theory still lacks some fea- 
tures which might be relevant for the dynamics at inter- 
mediate wavenumbers q ~ qcc- Now, we provide addi- 
tional evidence showing that this is the case. 

Figure 12(a) exhibits the simulation result for the tem- 
perature dependence of the a-relaxation times normal- 
ized by the one at q = q^ax, 'r,^ /"'"(^^x ' °^ density cor- 
relators (pg{t) for the density p = 2.83 molecules/nm^. 
MCT predicts the universal asymptotic power-law in- 
crease of the a-relaxation times [1]. This implies that 
the ratio of the a-relaxation times like the one shown in 
Fig. 12(a) should be temperature independent for T close 
to but above Tc, which is referred to as the a-relaxation- 
scalc coupling. However, Fig. 12(a) clearly demonstrates 
that this universal prediction of MCT is violated around 
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(and only around) the wavcnumber qcc- (For the other 
wavenumbers, it is seen that the a-relaxation-scale cou- 
pling holds fairly well, including T = 230 K which is 
below Tc « 234 K for the density considered.) Thus, 
the found temperature dependence of the ratio / t^^^ 
around q = qcc is beyond the implication of MCT, and 
the present theory, formulated within the framework of 
MCT, cannot be used to explain such an finding. Unfor- 
tunately, due to the presence of incoherent background, 
it is not clear whether experimental results for the coher- 
ent a-relaxation times reported in Ref. [12] also exhibit 
such a temperature dependence at q ^ qoc- (See also be- 
low for some other related simulation and experimental 
studies.) 

Combining the results shown in Figs. 4(a), 4(b) and 
12(a), one recognizes another interesting property of the 
unusual peak at g w gcc that it is more pronounced at 
lower density and at higher temperature. In view of the 
fact that the overall shape of the LW OTP molecule is 
well described as a sphere of the van dcr Waals radius 
rw = 0.37 nm (c/. Fig. 1), the dynamics at low-density, 
higli-tomperature regime is expected to be dominated by 
the spatial correlation of molecule's geometrical center, 
i.e., by S^'". This is in contrast to the high-density, 
low-temperature glassy regime where the dynamics is 
primarily determined by the cage effect manifested as 
the main peak in . Thus, one might conjecture that 



the unusual feature at (7 « goc is an inheritance from 
the low-density, high-temperature dynamics, and that 
this cannot be explained by the universal predictions of 
MCT since, up to the lowest temperature investigated, 
the dynamics at 5 « qoc might not have yet reached 
the asymptotic regime for which those MCT predictions 
are applicable. Such a possibility has been discussed in 
Refs. [29, 34], where the standard MCT scenario for the 
glass-transition dynamics was shown to be modified for 
some reorientational correlators due to precursor phe- 
nomena of a nearby type-A transition. 

The above conjecture, however, might not be appro- 
priate in the present case in view of the following re- 
sult. Figure 12(b) exhibits the temperature dependence 
of the plateau heights of the correlators (f>^{t). For 
T > Tc, the plateau height can be obtained from the von- 
Schweidler-law fit according to Eq. (6). Although Eq. (6) 
cannot be employed for fitting the correlators referring to 
T < Tc, we used it just to estimate the plateau heights 
for T = 230 K. (For T = 230 K, the correlators (f>^{t) 
exhibit well developed plateau regime as, e.g., shown in 
Fig. 3 of Ref. [23] , and the plateau heights can easily be 
estimated from such graphs. We confirmed that the so- 
estimated plateau heights are in good agreement with the 
ones based on the ad-hoc use of the von-Schweidler-law 
fit.) 

Since the plateau heights for T > Tc corresponds to the 
critical nonergodicity parameters /^'^ (c/. Eq. (6)), MCT 
predicts that they should be temperature independent. 
Furthermore, MCT predicts for T < Tc the increase of 
the plateau height, — f^'^ oc /i^ > 0. The results 
shown in Fig. 12(b) are, within the statistical errors, con- 
sistent with these universal predictions of MCT. Notice 
that this holds also for the intermediate wavenumbers 
Q ~ qcc- Thus, the plateau heights around q = qcc have 
already reached the asymptotic regime for which the uni- 
versal MCT description is adequate. Since the plateau 
heights also quantify the strengths of the a-relaxation 
processes, it is then difficult to imagine that only the a- 
rclaxation times have not yet reached the MCT asymp- 
totic regime, and the above conjecture introduced to ex- 
plain the finding in Fig. 12(a) might not be appropriate. 

The position q « gcc where the unusual feature we 
discussed occurs in LW OTP is compatible with the in- 
verse of its van der Waals radius {of. Sec. Ill A), i.e., it is 
connected to the overall size of the molecule. It is inter- 
esting to note that similar unusual peaks were found in a 
model for polymer around the intermediate wavcnumber 
q = 2TT/Rg, where Rg denotes the radius of gyration [53]: 
in this wavcnumber regime, a shoulder is discernible in 
the critical nonergodicity parameters and peaks are 
observable in the inverse of the critical amplitudes l/hg, 
the a-relaxation times Tg, and the stretching exponents 
(3q of the correlators which correspond to (p^{t) in the 
present paper. In particular, the ratio Tq/zg^^^ of the 
a-relaxation times at g' « 27r/i?g for this model exhibits 
the same temperature dependence as the one shown in 
Fig. 12(a). A similar temperature dependence of the ra- 
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tio Tq/Tq^,^^ at intermediate q range 0.4 (/max) was also 
found in the coherent neutron scattering results for a real 
polymer system [54]. Thus, further investigations are 
necessary for a comprehensive understanding of the un- 
usual features at intermediate wavenumbers as observed 
in simulation and experimental results for molecular and 
polymer systems. 
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